struct particula
{
    ///Posicion
    double x,y,z;
    ///Velocidad
    double vx,vy,vz;
    ///Diametro
    double diam;
    ///Masa
    double mass;
};

struct parametros {
    ///Lados de la caja
    double lx, ly, lz;
    ///Numero de particulas
    int npart;
    ///Paso de la simulacion
    int nacc;
    ///Tiempo
    double t;
    ///Intensidad del campo
    double u0;
    ///Volumen de la caja
    double vol;

};

void importBinary (ifstream *filein, particula *part,  parametros *param)
{
    double x,y,z,vx,vy,vz,diam,fx,fy,fz;
    double lx,ly,lz;
    double t;
    int npart,nacc;
    (*filein).read(reinterpret_cast<char *>(&t), sizeof(double));
    (*filein).read(reinterpret_cast<char *>(&nacc), sizeof(int));
    (*filein).read(reinterpret_cast<char *>(&npart), sizeof(int));
    (*filein).read(reinterpret_cast<char *>(&lx), sizeof(double));
    (*filein).read(reinterpret_cast<char *>(&ly), sizeof(double));
    (*filein).read(reinterpret_cast<char *>(&lz), sizeof(double));
    param->t=t;
    param->lx=lx;
    param->ly=ly;
    param->lz=lz;
    param->nacc=nacc;
    param->npart=npart;
    for (int j=0;j<npart;j++)
    {
        (*filein).read(reinterpret_cast<char *>(&x), sizeof(double));
        (*filein).read(reinterpret_cast<char *>(&y), sizeof(double));
        (*filein).read(reinterpret_cast<char *>(&z), sizeof(double));
        (*filein).read(reinterpret_cast<char *>(&vx), sizeof(double));
        (*filein).read(reinterpret_cast<char *>(&vy), sizeof(double));
        (*filein).read(reinterpret_cast<char *>(&vz), sizeof(double));
        (*filein).read(reinterpret_cast<char *>(&fx), sizeof(double));
        (*filein).read(reinterpret_cast<char *>(&fy), sizeof(double));
        (*filein).read(reinterpret_cast<char *>(&fz), sizeof(double));
        (*filein).read(reinterpret_cast<char *>(&diam), sizeof(double));
        part[j].x=x;
        part[j].y=y;
        part[j].z=z;

        part[j].vx=vx;
        part[j].vy=vy;
        part[j].vz=vz;
        part[j].diam=diam;
        part[j].mass=pow(part[j].diam,3.0);
    }

    param->vol=(param->lx)*(param->ly)*(param->lz);

}


double import(ifstream *MiArchivo, particula *part,  parametros *param)
{
    char cdummy[80];
    double t;
    *MiArchivo>>cdummy>>t>>cdummy>>param->nacc>>cdummy>>param->npart
        >>cdummy>>param->lx>>param->ly>>param->lz;
    param->t=t;
    double rr;
    ///Importamos las posiciones y velocidades
    for(int i=0;i<param->npart;i++) {
        ///Posiciones
        *MiArchivo>>rr;
        part[i].x=rr;
        *MiArchivo>>rr;
        part[i].y=rr;
        *MiArchivo>>rr;
        part[i].z=rr;

        ///Velocidades
        *MiArchivo>>rr;
        part[i].vx=rr;
        *MiArchivo>>rr;
        part[i].vy=rr;
        *MiArchivo>>rr;
        part[i].vz=rr;

        ///Fuerzas
        *MiArchivo>>rr;
        *MiArchivo>>rr;
        *MiArchivo>>rr;


        *MiArchivo>>rr;
        part[i].diam=rr;
        part[i].mass=pow(part[i].diam,3.0);
        param->vol=(param->lx)*(param->ly)*(param->lz);
        //box->boundary(i);
    }

    return t;
}

inline void klingemberg(double dist, double u0, double rcut, double s1, double s2, double *frep, double *epoten)
{
    double dummy;
    double smed=0.5*(s1+s2);
    if (dist<smed+rcut)
    {
        dummy=u0*exp(-100*(dist-smed));
        ///Para las fuerzas
        *frep=-3.0*dummy;
        ///Para el potencial
        *epoten=0.00375*dummy;
    }
    else
    {
        *frep=0.;
        *epoten=0.0;
    }
}

inline void RShiftedLJ(double dist, double sigma, double rcut, double s1, double s2, double *frep, double *epoten)
{
    double dummy;
    double smed=0.5*(s1+s2);
    double r0,roff,frac2,frac6;
    r0=smed-sigma;
    double rc=smed+2.5*sigma/r0;
    if (dist<rc)
    {

        roff=dist-r0;
        frac2 = (sigma/roff)*(sigma/roff);
        frac6 = frac2*frac2*frac2;

        ///Para las fuerzas
        *frep   = -48.0 * frac6*(frac6 - 0.5) / (roff);
        ///Para el potencial
        *epoten= 4.0*(frac6*frac6-frac6);
    }
    else
    {
        *frep=0.;
        *epoten=0.0;
    }
}

inline void magnetic(double dist, double dz, double u0, double *fpartxy, double *fpartz, double *epoten)
{
    double Co, Co2;
    double dist3= dist*dist*dist;
    double dist4= dist3*dist;
    double u30dis=3.0*u0/dist4;
    ///Calculamos la fuerza entre iii y jjj
    Co=dz/dist;
	Co2=Co*Co;
	*fpartxy= u30dis*(5.0*Co2-1.0);
	*fpartz= u30dis*(5.0*Co2-3.0);
	*epoten = u0*(1.0-3*Co2)/(dist3);
}


double interaccion(particula *part, parametros param, double *spxy, double *spxz, double *spyz, double *spzz)
{
	///Si calcStress = true, se calcula la parte correspondiente al potencial
	///del tensor de esfuerzos microscopicos y se actualiza en el objeto caja
	double distancia2(particula part1,particula part2, double *dx, double *dy, double *dz, parametros param);
	double massi,massj;
	double fMagxy, fMagz, fMagE;
	double frep, erep;
	double dist, dist2, ePoten;
	double fx,fy,fz;
	///Tensor de esfuerzos
	double sxy, sxz, syz, szz;
	///Radio de las particulas interactuantes
	double s1;
	double s2;
	///Campo magnetico
	double u0;
	double dx,dy,dz;
	double rcut=2.7;
	double sigma=0.1;

	///Inicializamos la energia potencial
	ePoten=0.;
	///Inicializamos el tensor de esfuerzos
	sxy=sxz=syz=szz=0.0;
	fMagxy=fMagz=fMagE=0.0;
	///Recorremos todos los pares de particulas
	for(int iii=0;iii<param.npart-1;iii++)
	{
		for(int jjj=iii+1;jjj<param.npart;jjj++)
		{
			dist2= distancia2(part[iii],part[jjj], &dx,&dy, &dz, param);
			dist=sqrt(dist2);
			s1=part[iii].diam;
			s2=part[jjj].diam;
			massi=part[iii].mass;
			massj=part[jjj].mass;
			u0=param.u0*massi*massj;

			magnetic(dist,dz, u0, &fMagxy, &fMagz, &fMagE);
			//klingemberg(dist, u0, rcut, s1, s2, &frep, &erep);
			RShiftedLJ(dist, sigma, rcut, s1, s2, &frep, &erep);

			///Calculamos la fuerza entre iii y jjj
			fx= (fMagxy+frep)*dx/dist;
            fy= (fMagxy+frep)*dy/dist;
            fz= (fMagz+frep)*dz/dist;
//            fx= (fMagxy)*dx/dist;
//            fy= (fMagxy)*dy/dist;
//            fz= (fMagz)*dz/dist;

            sxy+=dx*fy;
            sxz+=dx*fz;
            syz+=dy*fz;
            szz+=dz*fz;

            ///Acualizamos la energia potencial
			ePoten+=fMagE;
		}
	}
    ///Actualizamos la parte del potencial del tensor de esfuerzos microscopico
    *spxy=sxy/param.vol;
    *spxz=sxz/param.vol;
    *spyz=syz/param.vol;
    *spzz=szz/param.vol;


	ePoten/=param.npart;
	return ePoten;
}

double distancia2(particula part1,particula part2, double *dx, double *dy, double *dz, parametros param)
{
	///Calcula la distancia entre las particulas i y j
	///Consideramos condiciones de contorno periodicas en el plano xy
	///pero no en el z donde se encuentran las paredes

	double dis2=0;

    *dx=part2.x-part1.x;
	while(*dx>0.5*param.lx){*dx-=param.lx;}
    while(*dx<-0.5*param.lx) {*dx+=param.lx;}

    *dy=part2.y-part1.y;
    while(*dy>0.5*param.ly) {*dy-=param.ly;}
	while(*dy<-0.5*param.ly) {*dy+=param.ly;}

    *dz=part2.z-part1.z;
    #ifndef CONFINED
    while(*dz>0.5*param.lz) {*dz-=param.lz;}
	while(*dz<-0.5*param.lz) {*dz+=param.lz;}
    #endif

	dis2=(*dx)*(*dx)+(*dy)*(*dy)+(*dz)*(*dz);
	return	 dis2;
}

